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Abstract. The hydrodynamic stability of accretion disks is considered. The particular question is whether the combined action 
of a (stable) vertical density stratification and a (stable) radial differential rotation gives rise to a new instability for nonax- 
isymmetric modes of disturbances. The existence of such an instability is not suggested by the well-known Solberg-H0iland 
criterion. It is also not suggested by a local analysis for disturbances in general stratifications of entropy and angular momentum 
which is presented in our Section 2 confirming the results of the Solberg-H0iland criterion also for nonaxisymmetric modes 
within the frame of ideal hydrodynamics but only in the frame of a short-wave approximation for small m. As a necessary 
condition for stability we find that only conservative external forces are allowed to influence the stable disk. As magnetic forces 
are never conservative, linear disk instabilities should only exist in the magnetohydrodynamical regime which indeed contains 
the magnetorotational instability as a much-promising candidate. 

To overcome some of the used approximations in a numerical approach, the equations of the compressible adiabatic hydro- 
dynamics are integrated imposing initial nonaxisymmetric velocity perturbations with m = 1 to m = 200. Only solutions 
with decaying kinetic energy are found. The system always settles in a vertical equilibrium stratification according to pressure 
balance with the gravitational potential of the central object. 
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1. Introduction: The Solberg-Hoiland criterion 

Keplerian disks are, according to the Rayleigh criterion 

dj 2 

with j = R 2 Q being the angular momentum per unit mass 
and the distance R from the rotation axis, hydrodynamically 
stable configurations. They are unstable configurations, how- 
ever, under the influence of a (weak) magnetic field if the 
disk plasma has a sufficiently high conductivity (Balbus & 
Hawley 1991; Balbus 1995; Brandenburg et al. 1995; Ziegler & 
Rudiger 2000). The question has now been formulated whether 
also other, nonmagnetic influences exist which in combination 
with the basic (negative) shear dfi / dR allow instabilities in the 
linear regime. We remember three recent discussions. Urpin & 
Brandenburg (1998) consider the destabilizing action of any 
vertical dependence of Q on z which may more than compen- 
sate the stabilization of the positive radial gradient of j. 

There is also the suggestion of Klahr & Bodenheimer 
(2001) that the negative radial entropy stratification in thin 
Keplerian disks may act in the sense of a destabilization. Note 
that with the standard-alpha description of accretion disks one 



(2) 



finds the positive value 
dS_ _ 1 C„ 
dR ~ 2~R' 

indicating, however, stability. Finally Richard et al. (2001) and 
Brandenburg & Dintrans (2001) comment on the idea that the 
vertical gradient of the density which always dominates accre- 
tion disk stratifications in combination with positive dj /dR 
may destabilize the disk in the hydrodynamic regime. The ex- 
tensive numerical simulations in the local shearing box system 
have shown that Keplerian disks (without vertical gradient of 
the angular velocity) are completely stable to linear and non- 
linear instabilities (Balbus et al. 1996). 

Most of the mentioned discussions are covered by the 
Solberg-H0iland conditions for dynamic stability which, for 
the rotation studied here, take the form 



1 dj 2 1 
it' 3 dR C p p 



dp (dj 2 dS 
dz V dR dz 



dj 2 dS\ 



(3) 



(4) 
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(Tassoul 2000). If dp/dz < (as is usual for accretion disks) 
the latter relation turns into 

dj 2 dS dj 2 dS 



dR dz 



> 



dz dW 



(5) 
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which for Keplerian disks with Q = Q(R) oc R 15 simply 
leads to the usual Schwarzschild criterion 



dS_ 

dz 



> 



(6) 



for stability. However, if the angular velocity fl depends on 
the vertical coordinate z then Eq. (|J) allows instability in the 
standard case of positive dj/dR and dS/dz in case they are 
large enough. 

Equation (^|) becomes 



R 3 dR C v 



dS_ 

dz 



1 dp dS 



> 



(7) 



pC p dR dR 

which with (g) and for dS/dR > ensures stability. It was 
adopted that — g z = g = f2 2 z > 0. 

In the opposite case of a negative entropy gradient one at 
least needs 



dS 



dR 



> 



P C P Q 2 



(8) 



\dp/dR\ 

in order to violate Eq. (^|) in the most optimistic case of vertical 
isentropy. 

If, on the other hand, the standard stratifications of accre- 
tion disks are j — j(R) and S — S(z), with dj 2 /dR > then 
the Solberg-H0iland criterion reduces to 



1 



z dS 
C p dz 



f>0. 

dz 



(9) 



(in the upper hemisphere). Obviously dS/dz > is a suffi- 
cient condition for stability of Kepler disks after the Solberg- 
Hpiland criterion also for the combined action of vertical den- 
sity stratification and differential Kepler rotation. The same 
result have been derived in earlier papers by Livio & Shaviv 
(1977), Abramowicz et al. (1984) and Elstner et al. (1989). 
The local analysis for ring-like disturbances by Elstner et al. 
demonstrates in all details that in the framework of an ideal 
hydrodynamics, i.e. for vanishing Shakura-Sunyaev alpha, the 
traditional Schwarzschild criterion dS/dz > ensures stability 
for the Kepler disk. If the disk can be considered as isothermal 
in the vertical direction then it follows 



dS_ 

dz 



1Z d log p 



dz 



(10) 



which is always positive for the typical density stratification, 

1. e. the disk should always be stable. 

2. Combined stability conditions 

In this Section a local stability analysis of the equations of 
ideal hydrodynamics in cylindrical coordinates (R, <fi, z) is pre- 
sented. Sound waves and nonaxisymmetric disturbances are in 
particular allowed to exist. 

The three components of the momentum equation are 



du R 
dt 



ur 



duR U& du R 
dR R d<j> 
1 dp 
~ ~~p"dR 



du R 



R 



9r, 



du 6 dus us du A du A usur 

^r +UR mi + ii^ +u *^z- + — 

1 dp 

= ~~p~Rd4> +9 ^ 

du z du z us du z du z 

-dT + UR jR + li^ + u > 

The mass conservation reads 



dz 



Idp 

p dz 



9z- 



(12) 
(13) 



dp d , , pur 19, . d , . „ ,„ ^ 
It + dR {pUR) + -R + R66^ + dl^ = ° (U) 



and the energy equation is 

dS us dS 



dS_ 
~dt 



ur 



dR R 



as 

dz 



0. 



(15) 



As usual p is the density, p is the pressure, S is the specific en- 
tropy, g is the vector of any acceleration. The unperturbed state 
is assumed to be ur^q = u z _ = 0, usfi = Rf2(R, z), p a = 
Po(R,z),p Q = p (R,z), S = S (R,z),g R = g R (R,z), 
gs — 0, g z = g z (R,z). The only non-zero equations of the 
system (|TT|)-([T5[) are equations ( pT| ) and (|l3]) which take the 
form 



9R 



1 dp 
Po dR 



Rfl 2 



1 dp 
Pa dz ' 



(16) 



Finally, the ideal gas equation of state is used p = IZ/ppT, 
dS = C v d\og(pp^ 1 ) with C v as the specific heat at constant 
volume, and 7 = C p /C v = 5/3. 

2.1. Linearized equations 

The perturbations in the force are neglected. The perturbations 
ur,ua,,u z , pi, pi to the basic state are assumed small and the 
linearized system (JTTJ)-(JT^) is 



du R 



n 



dui 



dt dcj) 
dus 1 d 



1 dpi pi dpo 

2nu * + 7odR-7o^i =0 ' 



,du d 



an 

it— u z 

dz 



1 dpi 



0. 



p R 

du z , _ P±^Po _ Q 

dt dcj) po dz Pq dz 
1 dpi ^ du R ^ ur ^ du^ ^ 1 dus | 

Po dt 



dR R dz R dcj) 



Q dpi dlogpo 



Po d(j> 

dSi n dSi ds 
-dT + n ^ + UR mi 



dR 



-Ur 



dlogpo 

dz 



dS Q . 
u z -^— = 0. 

dz 



(17) 



(18) 



(19) 



u z =0, (20) 



(21) 



In the frame of our local approximation the equations are 
considered only in small volume near some reference point 
(i?o, 4>o, zq). The coefficients in the equations are thus constant 
taken at (i?o, <fio, zq). As usual both the short-wave approxima- 
tion and the small-m approximation, i.e. 



( U ) \k R R\^>m, \k z z\>m 



(22) 
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are applied (cf. Meinel 1983). In this approximation all 
perturbations can be expressed by the modal representation 
exp[i(fcfli? + m<j> + k z z — cut)] and the equations take the final 
form 



n\ on i -u Pl P 1 dpo n 

U) — mll)UR — lilUA, + IKr 2~n77 = 

Po Po v& 



u> ■ 



m(l)u<f, + ^-^{R 2 i2)u R + R^-u z 
mil)u z + ik z — - ^k-j^- = 0, 



0. 



Po Po dz 



pi dlogpo 

uj - mil) 1 — u R 

Pa oR 



uj — mil) 7 — 

Po Po 



dlogpo 

— » 1 

dz 

- ik R u R ^ 



1 

"c v 



dS, 



o 



OR 

Generally, it is 



Ur 







dz 



= 0. 



1 dS 8 (po 

= ^ l0S U 



C v dR dR 



(23) 
(24) 
(25) 



ik z u z = 0, (26) 



(27) 



1 dS d ( po . 

a^r = ^ l0g U )<28) 



The epicyclic frequency k and the adiabatic sound speed c a 
are 



K 2 



J_df 

i? 3 dR' 



7" 



Po 
/—■ 

PO 



(29) 



The determinant of the above homogeneous fifth order system 
must vanish and the resulting dispersion relation is 



The stability of axisymmetric and nonaxisymmetric pertur- 
bations is thus defined by the same stability criterion in the ap- 
proximation used. The first factor in ( j30| ) describes the waves 
with u = mil and will not be considered below. The stability 
is defined by the factor D. The roots of D are 



(u - mil) 2 =E± \/E 2 - F. 



(34) 



According to (32) the coefficient E is real. The flow is thus 



always unstable for negative E and complex F. Positive E 
and real F are, therefore, the necessary conditions for stabil- 
ity. According to ([33]), F is real if and only if 



1 / dp Q dp 
dz dR 



Pi 



dp dp 
dR dz 



R 



dQ 2 



0. 



Inserting (B5h in ( 16 ), we find 



dz 



dg z 
dR 



(35) 



(36) 



as the immediate consequence. Without forces this relation is 
always fulfilled. Any conservative force is a particular solution 
of (|36|). If - as it is in accretion disks - the gravity balances the 
pressure and the centrifugal force, then Eq. ( |3~5| ) is automati- 
cally fulfilled. Note that after the Poincare theorem for rotating 
media with potential force and fl = Q(R) both the density 
and the pressure can be written as functions of the generalized 
potential so that ( |3~5| ) is always fulfilled. 

Generally, the magnetic field is not conservative and can 
never fulfill the condition (|36|). This is the basic explanation 
for the existence of the magnetorotational instability which is 
driven by (weak) magnetic fields. 



(uj - mil) ■ D = 
with 

D = (uj- mil) 4 - 2E{uj - mil) 2 + F. 
Here we have introduced the expressions 



E = 

and 

F = 



(k R + k 2 )c^ c + 



1 dpo dpo 
p 2 dR dR 



1 dp dp 

Pq dz dz 



k R dp 
po dz 



kRcl c d 



log 



k z dp 
po dR 

7 

Po 



7 dz 
j_ ( dpodpo _ dp dp 
p 2 \ dz dR dR dz 
dfl 2 { , , 9 1 dpo dpo 

7 Q 



7 dR VPo 



(30) 



(31) 



(32) 



2 



k 2 r 2 



R ~d~z~ ' klikzcI 



1 dp dp Q 
Pq dz dz 
kR&Po _ k z dp 
po dz po dR 

(33) 



For axisymmetric modes we simply have to replace uj — mfl 
by uj. This procedure influences the frequencies rather than the 
stability. 



2.2. Sufficient conditions for stability 

Now we can suppose that the necessary condition ( [3~5| ) is ful- 
filled. 



The flow is stable if both roots in (34) are real and positive. 
The sufficient conditions for stability are therefore 



E > 0, 
Inserting 



< F < E 2 . 
into (M) we find F to be positive if 



(37) 



^ 2 

1,2 z 



where 



kR 2 dpo d 
k z 7Po dz dR 
1 dpo 

k z c a.cPo @ Z 



log 



dpo r,dil 2 dpo 



dz 



R 



dz dR 



N 2 = _J_dPod_ los fpo 

7Po dz dz VPo 



> 0,(38) 



(39) 



and 



1 dp d 



7Po dR dR 



Po 



Po 



(40) 



are components of the Brunt-Vaisala frequency. In our short- 
wave approximation (see (|2|)) we can neglect the last term in 
the left-hand side of (Q). The expression (§8|) is thus a simple 
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quadratic expression in k^Jk z . Two conditions will ensure its 
positivity: i) the expression is positive for some value kn/k z 
and ii) the expression has no real roots i.e. the discriminant is 
negative. The first condition is fulfilled if 



Nr + N 2 + k 2 > 0. 



(41) 



This is the Schwarzschild criterion for stability under the pres- 
ence of rotation. The second condition leads to 



d Po ( 2 d , ( Po 

M K ^ log U 



on 2 d fpo 



< 0. (42) 



Equations ( |41| ) and (|42| ) are exactly equivalent to the Solberg- 
H0iland conditions (||) and (Q). 

The Schwarzschild criterion ( |4l| ) ensures that E > and 
the short-wave approximation ensures that F < E 2 and these 
relations do not yield any additional conditions. 

The Solberg-H0iland conditions (|j) and (Q) can also be ob- 
tained in the Boussinesq approximation. Goldreich & Schubert 
(1967) have already formulated the dispersion relation within 
the frame of the Boussinesq approximation. Nevertheless, our 
more general consideration allows to find the new necessary 
condition for stability (|3T 

2.3. Limiting cases 

Let us consider two interesting limiting cases. Without rotation 
the Schwarzschild criterion (||) takes the form 

N 2 R + N 2 > 0, (43) 

where Nr and N z are given by © and ©. The stratification 
is thus always stable if po and So have opposite stratifications 
- which indeed is the standard case. Nevertheless, Eq. ( p3| ) al- 
lows the stability also for the case when po and So have the 
same stratifications in one direction and opposite stratifications 
in another. Note, that according to ( p3[ ) the situation with pres- 
sure stratification without any density stratification is always 
unstable. 

Without density stratification the sufficient stability condi- 
tion (pi) for a rotating flow is 



K 2 > 



1 



Po C ac 



dpo 
dR 



1 



Po C ac 



dpo 
dz 



(44) 



This is the classical Rayleigh criterion ([!]) for stability general- 
ized to the compressible case. Obviously, the right-hand side is 
positive and the criterion (Q) is more restrictive than the stan- 
dard one. 



finite disturbances. Along this way also finite-amplitude distur- 
bances can be considered within the frame of a nonideal hydro- 
dynamics. 

The above analytical study has shown that stratified disks 
are stable; under the assumption of a small scale in perturba- 
tions compared with the disk dimensions. We present a few nu- 
merical integrations of the hydrodynamic equations for the un- 
restricted case. The computations are based on the work of Arlt 
& Riidiger (2001) which required a stable, stratified Keplerian 
disk as a prerequisite. We adapt these simulations with adia- 
batic evolution and relatively strong nonaxisymmetric pertur- 
bations. 

The setup for our three-dimensional simulations applies a 
global, cylindrical computational domain with full azimuthal 
range, dimensionless radii from 4 to 6, and a vertical extent of 
2 density scale heights on average, which is —1 to +1 in di- 
mensionless units. The ZEUS-3D code (developed by Stone & 
Norman 1992a,b; Stone et al. 1992) was used for the integra- 
tion of the hydrodynamics. The modifications to the code are 
very close to those given in Arlt & Riidiger (2001); the induc- 
tion equation is naturally not integrated here. 



dp 

at 



+ div(pit) = 0, 



du 

P ~dt + P( uV > u = + p9 ' 



de 



+ div(ett) = — pdivu, 



(45) 



(46) 



(47) 



where it, p, e, and p have the usual meanings of velocity, den- 
sity, energy density, and pressure. The source of gravitation is 
that of a point mass Al = 10 5 at R = 0. The simulation does 
not include self-gravity in the disk. The polytrophic exponent is 
7 = 5/3. The original ZEUS code provides artificial viscosities 
for improved shock evolution; we have put these terms to zero 
here. The advection interpolation is the second order van-Leer 
scheme. 

The conditions for the vertical boundaries are stress-free; 
no matter can exit the computational domain in vertical direc- 
tion and the vertical derivatives of ur and u$ vanish. The radial 
boundaries are also closed for the flow, and the radial deriva- 
tive of u z vanishes. For the azimuthal flow, we have to adopt 
a modified boundary condition though. A Keplerian rotation 
profile is maintained into the boundary zones, based on the last 
zone of the integration domain. If the first azimuthal velocity 
of the computational domain is denoted by ul\ we use 



3. Numerical simulations 

The main shortcoming in our presentation is the short-wave 
assumption ( p2[ ) under which all the above considerations are 
only valid. Without this assumption, Fourier modes are no 
longer solutions of the linearized equations and we have to 
switch to much more complicated mathematics involving in- 
tegral equations (see Riidiger & Kitchatinov 2000). We prefer 
to present numerical simulations for the stability of (isother- 
mal) density-stratified Keplerian disks under the influence of 



(48) 
(49) 



The outer radial boundary is treated accordingly for Ui B+1 ' and 
^(n+2) rp ne az j mutna i boundary conditions are periodic. 

A rough representation of a stratified disk with radius- 
independent density scale-height, was used for the initial con- 
ditions. We leave this configuration for free development under 
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Table 1. Overview of simulation configurations. The quantity 
m denotes the azimuthal mode of the initial perturbation. 



Model 


Resolution 


m 


I 


31 x 61 x 351 


1 


II 


71 x 71 x 351 


1 


III 


31 x 61 x 701 


10 


IV 


31 x 61 x 701 


200 



u z = Asinmip, 
uji = A cos(m + 1) 



1 



7(7 - 1) 



where the sound speed c ac is 7% of the Keplerian velocity Uk, 
and the polytrope exponent is 7 = 5/3. 

Figure 1 shows the evolution of the directional kinetic en- 
ergies for the z- and r-components. A gradual decay of fluid 
motion is observed. Velocity fluctuations can also be measured 
in terms of the correlation tensor element Qr^ which is non- 
dimensionalized with the average pressure giving the Shakura- 
Sunyaev parameter ass (cf. Arlt & Riidiger 2001). Figure 2 
shows the evolution of ass during 30 rotation periods at the 
inner disk radius. A clear relaxation of the flow is seen. The 
initially strong velocity fluctuations settle the correct density 
stratification for the gravitational potential of the point mass 
(central star). The fluctuations in vertical direction reach ±14% 
of the Keplerian velocity Uk- Those in the radial direction even 
reach —27%. They are thus not small at all, and even a nonlin- 
ear instability might be noticed if existing. 

We can read an approximate decay time td from Fig. 1 . On 
average, it amounts to about 20 orbital periods or td = 3.18 
in non-dimensional units. As the main relaxation motions are 
vertical, we use the density scale-height for the length scale and 
obtain a viscosity of 0.03 caused by the finite resolution of the 
numerical grid. In terms if Reynolds numbers, we achieve 



Re 



H 



1520 



10 5 t 



the influence of the gravitational potential. The initial condi- 
tions also involve a non-axisymmetric velocity perturbation of 
the form 



j 1 [i kinetic in z-direction 



(50) 
(51) 




The amplitude A is - in terms of the Keplerian velocity in the 
middle of the annulus, Uk - about 7 • 10 -4 Z7k. We have run 
two models with m = 1 and different resolutions and an ad- 
ditional integration with rn = 10 in order to test short-wave 
perturbations. The models are summarized in Table |l} The ini- 
tial configuration is isothermal with an energy density of 



10 20 30 40 

TIME (ORBITS) 



50 



Fig. 1. Kinetic energies in the vertical and radial components of 
the flow from the run with 31 x 61 x 351 grid points (Model I). 
The solid line is the energy of the redirection, the dashed line 
that of the zd-direction. Times are given in revolutions of the 
inner boundary of the computational domain 



(52) 



5x10" 




(53) 



near the inner edge of the annulus where the Keplerian velocity 
is U la = 158. The density-scale height at the inner edge is 
Hp = 0.33 in the initial configuration. 

A second run with significantly higher resolution - in ver- 
tical direction in particular - is shown in Fig. 3. No change in 
the decaying nature of the perturbation is found. 

An increase in azimuthal wave number of the perturba- 
tion does not provide instability either. Figure [I] is the result 



-5x10" 



20 30 40 

TIME (ORBITS) 

Fig. 2. Shakura-Sunyaev parameter ass from the same model 
as in Fig. 1 (Model I). Times are again in orbital periods of the 
inner boundary of the domain 

of a simulation with a perturbation of very high wave number, 
u r = A sin 2000, u z — A cos 2010, and doubled azimuthal 
resolution. The gas again settles towards an equilibrium state 
and very low kinetic energies in the radial and vertical compo- 
nents. 

4. Conclusions 

We have shown in Section 1 that the Solberg-H0iland condi- 
tions (fy and (0) answer the question whether the combined 
action of stable (radial) shear and stable (vertical) density strat- 
ification can produce a new instability. The clear answer is No, 
but uncertainties remain about validity of the Solberg-H0iland 
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jiViiVi't kinetic in s " directlon 



10 z 




kinetic in r-direction 




10 15 20 
TIME (ORBITS) 



25 



Fig. 3. Kinetic energies in the vertical (dashed) and radial 
(solid) components of the flow. The resolution of the model 
is 71 x 71 x 351 grid points (Model II). Times are again in 
orbital periods of the inner boundary of the domain 



10 



10 z 



i'ljiii,','', kinetic in z— direction : 

-'lnl'l'll'il'" ii 1 ' 
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kinetic in r-direction 




10 15 20 25 30 
TIME (ORBITS) 



Fig. 4. Kinetic energies in the vertical (dashed) and radial 
(solid) components for the model with m = 200 in the initial 
perturbation (Model IV). The resolution of the computational 
domain is 31 x 61 x 701 



criterion for the case of nonaxisymmetric disturbances of the 
compressible accretion disk material. Indeed, in Section 2 we 
were able to reproduce the Solberg-H0iland criterion in the 
well-known formulation but only with a short-wave approx- 
imation and a small-m approximation applied formulated in 

Sound waves have been included into the consideration. 
Our dispersion relation ( pl| ) is of fourth order contrary to the 
dispersion relation of second order resulting in the Boussinesq 
approximation (Goldreich & Schubert 1967) and contrary to 
the dispersion relation of fifth order resulting in nonideal hy- 
drodynamics (Abramowicz et al. 1984, 1990). Along this way a 
new necessary condition for stability results which requires the 



external force which balances pressure and centrifugal force to 
be a conservative one, i.e. to possess a potential (cf. (j35|)). 

It is interesting to note that in the Boussinesq approx- 
imation for nonideal fluids the Solberg-H0iland criterion 
changes to the Goldreich-Schubert-Fricke criterion (Goldreich 
& Schubert 1967; Smith & Fricke 1975; Urpin & Brandenburg 
1998) which no longer allows the existence of a vertical shear 
for stability, i.e. dfl jdz = 0. And, indeed, even in thin accre- 
tion disks there is always such a (weak) vertical shear dfi/dz 
the meaning of which remains open. We have thus numerically 
simulated the stability of such accretion disks using the ZEUS- 
3D code. The numerical integration of the nonlinear adiabatic 
equations of the fully compressible hydrodynamics with ini- 
tially nonaxisymmetric velocity perturbations does not lead to 
the onset of instability although the disk, of course, does pos- 
sess a small vertical shear. An overview about the simulated 
configurations can be found in Table 1 which shows that we 
indeed have also considered very high values of the azimuthal 
"quantum" number m which do not fulfill the short-wave con- 
dition (H). 

Based on our numerical simulations, we thus cannot sup- 
port findings about a fast, hydrodynamic instability in strati- 
fied disk flows. Similar simulations including magnetic fields 
showed the magneto-rotational instability with growth rates 
near O.bP^. Our hydrodynamic runs cover 30-50 orbits of 
the inner disk and show only decay in velocity fluctuations. 
Angular momentum transport hovers at ass ~ ±10~ 6 ; the tem- 
poral average in Fig. 2 is 5 • 10 ~ 8 . 
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